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The accumulation of deleterious mutations is driven by rare fluctuations which lead to the loss of 
all mutation free individuals, a process known as Muller's ratchet. Even though Muller's ratchet 
is a paradigmatic process in population genetics, a quantitative understanding of its rate is still 
lacking. The difficulty lies in the nontrivial nature of fluctuations in the fitness distribution which 
control the rate of extinction of the fittest genotype. We address this problem using the simple 
but classic model of mutation selection balance with deleterious mutations all having the same 
effect on fitness. We show analytically how fluctuations among the fittest individuals propagate to 
individuals of lower fitness and have a dramatically amplified effects on the bulk of the population 
at a later time. If a reduction in the size of the fittest class reduces the mean fitness only after 
a delay, selection opposing this reduction is also delayed. This delayed restoring force speeds up 
Muller's ratchet. We show how the delayed response can be accounted for using a path integral 
formulation of the stochastic dynamics and provide an expression for the rate of the ratchet that 
is accurate across a broad range of parameters. 



By weeding out deleterious mutations, purifying selection acts to preserve a functional genome. In sufficiently 
small populations, however, weakly deleterious mutations can by chance fix. This phenomenon, termed Muller's 
ratchet (Felsenstein 1974 Muller |1964 ), is especially important in the absence of recombination and is thought 
to account for the degeneration of Y-chromosomes (Rice} 1987) and for the absence of long lived asexual lineages 



(Lynch et al. 



1993). 



A click of Muller's ratchet refers to the loss of the class of individuals with the smallest number of deleterious 
mutations. To understand the processes responsible for such a click, it is useful to consider a simple model of 
accumulation of deleterious mutations with identical effect sizes illustrated in Fig. [T] Because of mutations, the 
population spreads out along the fitness axis, which in this model is equivalent to the number of deleterious mutations 
in a genome. The population can hence be grouped into discrete classes each characterized by the number of deleterious 
mutations. Mutation carries individuals from classes with fewer to classes with more mutations, hence shifting the 
population to the left. This tendency is opposed by selection, which amplifies fit individuals on the right, while 
decreasing the number of unfit individuals on the left. These opposing trends lead to a steady balance, at least in 
sufficiently large populations. However, in addition to selection and mutation, the distribution of individuals among 
fitness classes is affected by fluctuations in the number of offspring produced by individuals of different classes, i.e., 
by genetic drift. Such fluctuations are stronger (in relative terms) in smaller populations and in particular in classes 
which carry only a small number of individuals. When mutation rate is high and selection is weak, the class of 
individuals with the smallest number of mutations (k = in Fig. [IJ contains only few individuals and is therefore 
susceptible to accidental extinction - an event that corresponds to the "click" of Muller's ratchet. 

Despite the simplicity of the classic model described above, understanding the rate of the ratchet has been a 



challenge and remains incomplete (Etheridge et al. 2007 Gessler 1995 GORDO and Charlesworth| |2000 
Higgs and Woodcock[ |1995[ |Jain[ |2008| |Stephan et ^ |1993| |Stephan and Kim[ |2002[|Waxman and Loewe 



2010). Here, we revisit this problem starting with the systematic analysis of fluctuations in the distribution of the 



population among different fitness classes. We show that fitness classes do not fluctuate independently. Instead, there 
are collective modes affecting the entire distribution which relax on different time scales. Having identified these 
modes, we calculate the fluctuations of the number of individuals in the fittest class and show how these fluctuations 
affect the mean fitness. Fluctuations in mean fitness feed back on the fittest class with a delay and thereby control 
the probability of extinction. These insights allow us to arrive at a better approximation to the rate of the ratchet. In 



particular, we show that the parameter introduced in the earlier work (Gordo and Charlesworth, 2000 Haigh 



1978 



Stephan et al. 1993 ) to parameterize the effective strength of selection in the least loaded class is not a constant 



but depends on the ratio of the mutation rate an d the effect size of mutations. We use the path integral representation 
of stochastic processes, borrowed from physics (Feynman and HiBBS 1965) to describe the dynamics of the fittest 
class and arrive at an approximation of the rate of Muller's ratchet that is accurate across a large parameter range. 
Understanding the rate of the ratchet is important, for example to estimate the number of beneficial mutations 
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required to halt the ratchet and prevent the mutational melt-down of a population (GOYAL et al. 2012 Lynch 



et al. 



1993 Pfaffelhuber et al. 



mutations we refer the reader to Charlesworth 



2012[) (for an in-depth and up-to-date discussion of the importance of deleterious 
,( [20T2| ). Furthe rmore, fluctuations of fitness distributions are a 
general phenomenon with profound implications for the dynamics of adaptation and genetic diversity of populations. 
Below we shall place our approach into the context of the recent studies of the dynamics of adaptation in populations 



with extensive non-neutral genetic diversity (Desai and Fisher 2007 Neher et al. 2010 Rouzine et al. 2003 



Tsimring et al. 1996). The study of fluctuations in the approximately stationary state of mutation selection balance 



which we present here is a step towards more general quantitative theory of fitness fluctuations in adapting populations. 



I. MODEL AND METHODS 



We assume that mutations happen at rate u and that each mutation reduces the growth rate of the genotype by 
s <C 1. Within this model, proposed and formalized by |Haigh individuals can be categorized by the number of 



deleterious mutations they carry. The equation describing the fitness distribution in the population, i.e., what part 
n k of the population carries k deleterious mutations, is given by 

^-n k = s(k - k)n k - un k + ran + ^fnir\ k (1) 

where k = N^ 1 ^2 k kn k i^2 k n k = N) and the last term accounts for fluctuations due to finite populations, i.e., 
genetic drift, and has the properties of uncorrelated Gaussian white noise with {r] k (t)rji{t')) — 5 k i8(t — t'). In the 
infinite population limit, this equation has the well known steady state solution n k = Ne~ x X k /k\, whe re A = u/s. A 
time dependent analytic solution of the deterministic model has been described in (Etheridge et al. |2007 ). 



Note that we have deviated slightly from the standard model, which assumes that genetic drift amounts to a 
binomial resampling of the distribution with the current frequencies N~ 1 n k . This choice would result in off-diagonal 
correlations between noise terms that stem from the constraint that the total population size is strictly constant. This 
exact population size constraint is an arbitrary model choice which we have relaxed to simplify the algebra. Instead, 
we control the population size by a soft constraint which keeps the population constant on average but allows small 
fluctuations of N. The implementation of this constraint is described explicitly below. We confirmed the equivalence 
of the two models by simulation. 



Computer simulations 

We implemented the model as a computer simulation with discrete generations, where each generation is produced 
by a weighted resampling of the previous generation. Specifically, 

^^^-^(l-s^V-^n (2) 




FIG. 1 Deleterious mutation selection balance. The population is distributed among classes of individuals carrying k deleterious 
mutations. Classes with few mutations grow due to selection (red arrows), but loose individuals through mutations (green 
arrows), while classes with many mutations are selected against but replenished by mutations. 
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Symbol 


Description 


N, u, s 


population size, mutation rate, and mutation effect 


n k {t), n k 


number of individuals in class k at time t, steady state value 


A = u/s, T = ts 


rescaled mutation rate and time 




population frequency in class k at time r, steady state value 


z T , z 


abbreviation for xq(t) and xo 


k 


mean fitness: k = kxt 


5xk, Sk 


deviations from steady state: Sxt = Xk — x~k, Sk = A — k 


Sx({z P }) 


path integral action depending on the path {z p } with < p < r. 


S*x(Zt,Z ), Z* 


extremal action and the associated path depending on the endpoints z t ,zq 


Sl(z) 


long time limit of Sx(z T , zo) with z = z T 


P T (Z T \Z()) 


propagator from zo to z T in time r 


P(z) 


steady state distribution of z 


7 


rate of the ratchet in units of s 


a 2 


variance of xo depending on Ns and A 


c 2 


rescaled variance of xo depending on A only: £ 2 = Nse x a 2 




k-ih component of right and left eigenvectors with eigenvalue — i 




projection of 5xk on tffj^ 


a 


parameter of the effective potential confining xo (traditionally a = 0.5-0.6) 



TABLE I List of symbols 



where W is the mean fitness W = CAT 1 £ fc (l - s) k n k and C = exp^A^ 1 - 1) is an adjustment made to the 
overall growth rate to keep the population size approximately at Nq. This specific discretization is chosen because it 

The simulation was 



1978) 



has exactly the same stationary solution as the continuous time version above (Haigh 
implemented in Python using the scientific computing environment SciPy ( |Oliphant |2007 ). If the parameter of the 
Poisson distribution was larger than 10 4 , a Gaussian approximation to the Poisson distribution was used to avoid 
integer overflows. 

To determine the ratchet rate, the population was initialized with its steady state expectation n^, allowed to 
equilibrate for 10 4 generations, and then run for further T = 10 8 generations. Over these 10 8 generations, the number 
of clicks of the ratchet where recorded and the rate estimated as clicks per generation. 

The source code of the programs, along with a short documentation, is available as supplementary material. In 
addition, we also provide some of the raw data and the analysis scripts producing the figures as they appear in the 
manuscript. 



Numerical determination of the most likely path. 

The central quantities in our path integral formulation of the rate of Muller's ratchet are i) the most likely path to 
extinction and ii) the associated minimal action S* x . To determine the most likely path to extinction of the fittest class, 
we discretize the trajectory into m equidistant time points pi between and r, where xo(0) = xq and Xo{p m ) = 0. For 
a given set of xo(pi), a continuous path xo(p) is generated by linear interpolation. For a given trajectory xq(p), we 
determine the mean fitness by solving the deterministic equations for Xk(p), k > 1. Note that in this scheme, the only 
independent variable is the path xo(p), all other degrees of freedom are slaved to Xo(p). From xo(p) and the resulting 



k(p), we calculate the action S\({xo(p)}) as defined in Eq. (28). S\({xo(p)}) is then minimized by changing the 
values of xo(pi), < i < m, using the simplex minimization algorithm implemented in SciPy ( |Oliphant |2007| ). To 



speed up convergence, the minimization is first done with a small number of pivot points (to = 4), which is increased 
in steps of 2 to m — 24. The total time r = 20 (in units of s _1 ) was used which is sufficiently large to make the result 
independent of r. The code used for the minimization is provided as supplementary material. 



II. RESULTS AND DISCUSSION 



Fluctuations of the size no of the least loaded class can lead to its extinction. In the absence of beneficial mutations 
this class is lost forever ( Muller 1964 ) , and the resulting accumulation of deleterious mutations could have dramatic 
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evolutionary consequences. Considerable effort has been devoted in understanding this process, and it has been noted 
that the rate at which the fittest class is lost depends strongly on the average number of individuals in the top class 
fio (Haigh 1978). Later studies have shown that the rate is exponentially small in n^s if fi s 3> 1 (Jain 2008). If 



tIqs is small, the ratchet clicks frequently and a traveling wave approach is more appropriate ( Rouzine et al.\ 2008). 
However, a quantitative understanding of the n^s 3> 1 regime is still lacking. 

Here, we present a systematic analysis of the problem by first analyzing how selection stabilizes the population 
against the destabilizing influences of mutation and genetic drift, and later use this insight to derive an approximation 
to the rate of Muller's ratchet. Before analyzing Eq. ([!]), it is useful to realize that it implies a common unit of 
time for the time derivative, the mutation rate and the selection coefficient which is of our choosing (days, months, 
generations, . . . ). We can use this freedom to simplify the equation and reveal what the important parameters are 
that govern the behavior of the equation. In this case, it is useful to use s _1 as the unit of time and work with the 
rescaled time r = ts. Furthermore, we will formulate the problem in terms of frequencies x k = N^ 1 n k rather than 
numbers of individuals, and obtain 



d 

Jr 



\)x k + \x k _ 1 + x j — % 



(3) 



where A — u/s is the dimensionless ratio of mutation rate and selection strength. In other words, A is the average 
number of mutations that happen over a time s" 1 (our unit of time). Note that A uniquely specifies the deterministic 
part of this equation and its steady state solution Xf. = \ k e~ x jk\. The stochastic forces are proportional to 1/V Ns. 
Again, the parameter combination (Ns) -1 has a simple interpretation as the variance of the stochastic effects accumu- 
lated over time s _1 . Other than through a prefactor determining the unit of time, any quantity governed by Eq. Ety 
can depend only on A and Ns. Hence it is immediately obvious that the ratchet rate cannot depend on uq = Ne 



alone, but has to depend on tiqs instead (Jain 2008). All times and rates in rescaled time units are denoted by greek 



letters, while we use roman letters for times and rates in units of generations. 

Before turning to the ratchet rate, we shall analyze in greater detail the interplay of deterministic and stochastic 
forces in Eq. A full time dependent analytic solution of the deterministic model was found in (Etheridge et al. 



2007). Below, we will present an analytic characterization of the stochastic properties of the system in a limit where 



stochastic perturbations are small. 




FIG. 2 Panel A: the covariance of the size of the fittest class 20 (0) with xo(t) a time r later. The normalized auto-correlation 
of xo increases with A. Panel B: the covariance of xo(0) with the mean fitness at time r in the past or future. One observes 
a pronounced asymmetry, showing that fluctuations of the fittest class propagate towards the bulk of the fitness distribution 
and results in delayed fluctuations. Simulation results are shown as dashed lines, theory curves are solid. In all cases, s = 0.01 
and Nxqs = 100. Note that time is measured in units of 1/s, which is the natural time scale of the dynamics. 



Linear stability analysis. 

In the limit of large populations, the fluctuations of x k around the deterministic steady state x k can be analyzed in 
linear perturbation theory. In other words, we express deviations from the steady state as Sx k = x k — x k and expand 
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the deterministic part of Eq. ^ to order Sxf, . This expansion 

d -Sx k = -k8x k + \5x k -i + x fe ^ (m - \)5x rn = ^ L km Sx m (4) 



dr 

m—0 rn 

defines a linear operator L km . A quick calculation shows that L km has eigenvalues «j = — i with i = 0, 1, 2 . . .. The 
right eigenvector corresponding to kq — is simply ip^ — x k , while the right eigenvectors for i > are given by 

ip\ % ) = x k -i - x k (5) 

where k numbers the coordinate of the vector. This is readily verified by direct substitution (Note that x,i = for 
i < 0). 

The eigenvector corresponds to population size fluctuations which in our implementation are a controlled by 
a carrying capacity. The eigenvalue associated with this mode in the computer simulation is large and negative and 
need not be considered here, see Methods. All other eigenvalues are negative, which is to say that x k is a stable 
solution. 

The eigenvectors for i > have an intuitive interpretation: Eigenvector i corresponds to a shift of a fraction of the 
population by i fitness classes downward. Since such a shift reduces mean fitness, the fittest classes start growing, 
and undo the shift. More generally, any small perturbation of the population distribution can be expanded into 
eigenvectors 8x k (r) = J2j ^ k ^ a j{ T ) an d the associated amplitudes %*(r) will decay exponentially in time with rate j 
(remember that the unit of time is s _1 ). Since the amplitudes are projections of 8x k onto the left eigenvectors of L mk , 
we need to know those as well. For kq — 0, the left eigenvector is simply (f>^ = 1, while the other left eigenvectors 
are given by 

(_-[\k-i p X\i-k 

and 0l = for k > i. With the left and right eigenvectors and the eigenvalue spectrum of the deterministic system 
on hand, we will now re- instantiate the stochastic part of the dynamics. 

^_5x k = ^2 L k m Sx m + Jjj^Vk (7) 

m 

Note that we approximated the strength of noise by its value at equilibrium. This approximation is justified as long 
as we consider only small deviations from the equilibrium. The full x k dependent noise term will be reintroduced 

later when we turn to Muller's ratchet. Substituting the representation of Sx k (r) — J2i a i( r ) an d projecting onto 

(i) 

the left eigenvector we obtain the stochastic equations for the amplitudes 

±a J (r) = -ja J (r) + J2^ ) Mv k (r) . (8) 

k 

Each noise term r\ k contributes to every ctj and induces correlations between the <Xj, but each amplitude can be 
integrated explicitely 

«i(r) = f dr'e-^-r ') ^\/1>(t') • (9) 

J— oo k V JVS 

The covariances of different amplitudes are evaluated in the Supplementary Information and found to be 

(a i (r)a j (r + Ar)) = e -—r E (10) 

J k 

However, we are not primarily interested in the covariance properties of the amplitudes of eigenvectors, but expect 
that the fluctuations of the fittest class and fluctuations of the mean fitness are important for the rate of Muller's 
ratchet and other properties of the dynamics of the population. To this end we express 6xq(t) and k(r) as 

Sx (r) - ^^« 3 (r) = -e-^ % (r) (11) 

j>0 j>0 

Sk(r) = - £ k4 j) aj (r) = -X>jto • (12) 

j>a,k j>o 



G 



Together with Eq. (10), we can now calculate the desired quantities. The calculations required to break down the 
multiple sums to interpretable expressions are lengthy, but straight-forward and detailed in the supplement. Below, 
we will present and discuss the results obtained in the supplement. 



Fluctuations of xo and the mean fitness. 

For the variance of the fittest class (Sxq) and more generally its auto-correlation, we find 

(fao(0)^(r)> = ^^^,r) (i3) 
G\(9, t) = e ™ 2 e--A(i+e-ne „ e -xe _ e -xoe-T + 1 

The variance of the fittest class (t = in the above expression) is therefore a 2 = ^C 2 (^) where C 2 (A) = Jq ^fG\(9, 0) 
is the standardized variance of the top bin, which depends only on A. For small A, it simplifies to ( 2 (\) ~ |A + C(A 2 ). 
This limit corresponds to xq close to 1 with only a small fraction of the population carrying deleterious mutations 
ii«A = u/s. The opposite limit of large A corresponds to a broad fitness distribution where the top class represents 
only a very small fraction of the entire population. In this limit, the leading behavior of the variance a 2 is ~ gfliog A , 
The full auto-correlation function is shown in Fig. [2j\ for different values of A and compared to simulation results, 
which agree within measurement error. In our rescaled units, the correlation functions decay over a time of order 
1, corresponding to a time of order 1/s in real time. More precisely, the decay time (in scaled units) increases with 
increasing A as log A. 

In a similar manner, we can calculate the auto-correlation of the mean fitness 

(5k(0)5k(r)) = ^ [ d0I x (0,r) 

Ns Jo (14) 

/ A (0, T ) = e -T e xe*e-T-x(i+e-ne (0 + xe (0g-T _ i) _ i)) 

which asymptotes to (ANsxq)^ 1 for large A at r = 0. It is hence inversely proportional to the size of the fittest class 
xq. For large A, xq represents only a tiny fraction of the population and fluctuations of the mean can be substantial 
even for very large N. This emphasizes the importance of fluctuations of the size of the fittest class for properties of 
the distribution. 

If fluctuations of the mean fitness k are driven by fluctuations of the fittest class Xq, we expect a strong correlation 



between those fluctuations (Etheridge et al. 2007). Furthermore, fluctuations of xq should precede fluctuations of 



the mean. These expectations are confirmed by the analytic result 



Ns 

where 



(Sx o (0)Sk(r)) = A j d6H x (6, r) (15) 



r ((9 _ 1)e _ r+e -^_A(l +e -) fl + e -r- e -^ r>0 
X{0 ' T >~\ - l) e e T Afl»-A(l+e^)fl + e -X6 T < g ^ 

This expression is shown in Fig. [2^3 for different values of A. The cross correlation (5xo(G)5k(r)) is asymmetric in 
time: With larger A, the peak of the correlation function moves slowly (logarithmically) to larger delays. This result 
is intuitive, since we expect that fluctuations in the fittest class will propagate to less and less fit classes and that the 
dynamics of the entire distribution is, at least partly, slaved to the dynamics of the top class. 

In all of these three cases, the magnitude of the fluctuations is governed by the parameter Ns, while the shape of 
the correlation function depends on the parameter A. Only the unit in which time is measured has to be compared 
to the strength of selection directly. 



The rate of Muller's Ratchet. 

The ratchet clicks when the size of the fittest class hits 0, and the rate of the ratchet is given by the inverse of the 
mean time between successive clicks of the ratchet. Depending on the average size Xq of the fittest class, the model 
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FIG. 3 An example of a click of the ratchet with TV = 5 x 10 , s — 0.01 and A = 10, corresponding to an average size of the 
fittest class n ~ 2269. Panel A: The distribution of xo averaged over the time prior to extinction. Panel B&C: The trajectory 
of xo(t), with the part of the trajectory that ultimately leads to extinction magnified in panel C. The final run towards Xo = 
takes a few time units, as expected from the results for the correlation functions which suggest a (rescaled) correlation time of 
~ log A. Note this time corresponds to a few hundred generations since s — 0.01. 



displays very different behavior. If Nsxq is comparable to or smaller than 1, the ratchet clicks often without settling to 
a quasi-equilibrium in between clicks. This limit has been studied in (Rouzine et al. 20081. Conversely, if Nsxq 3> 1 



ratchet clicks are rare and the system stays a long time close to its quasi-equilibrium state Xk- Such a scenario, taken 
from simulations, is illustrated in Fig. [3j Panel A shows the distribution of x$ prior to the click, while Fig. [3)3 shows 
the realized trajectory which ends at x = 0. Prior to extinction, Xq(t) fluctuates around its equilibrium value and 
large excursions are rare and short. The final fluctuation which results in the click of the ratchet is zoomed in on in 
Fig. [3p. Compared to the time the trajectory spends near xq, the final large excursion away from the steady is short 
and happens in a few units of rescaled time. Translated back to generations, the final excursion took a few hundred 
generations (s = 0.01 in this example). 

In rescaled time, the equation governing the frequency of the top class is 



dr 



Xq(t) = Sk(T)x (T) + 



go(r) 
Ns 



(17) 



where 5k(r) = X — fc(r). The restoring force 5k depends on x , as well as on the size of the other classes Xf.. For 
sufficiently large A, xq is much smaller than Xk with k > 1, such that the stochastic force is most important for x$. 
The dynamics of Xk, k > 1, is approximately slaved to the stochastic trajectory of xq(t). We can therefore try to 
find an approximation of 5k(r) in terms of Xq(t) only. The linear stability analysis of the mutation selection balance 
has taught us that the restoring force exerted by the mean fitness on fluctuations in xq is delayed with the delay 
increasing with oc log A. The latter observation implies that the restoring force on xq will mainly depend on the 
values of xo some time of order log A in the past. Such history dependence complicates the analysis, and this delay 
has been ignored in prev ious analysis, which assumed that 5k(r) depends on the instantaneous value of xq(t) via 



5k(r) = a(l - x (r)/a;o) (Gordo and Charlesworth 2000 Jain 2008 Stephan et al 



was chosen ad hoc between 0.5 and 0.6. This restoring force is akin to an harmonic potentia" 



1993). The parameter a 



centered around xq and 



the stochastic dynamics is equivalently described by a diffusion equation for the probability distribution P(xq,t). 



± P(z T ) = — 

dr 1 ' ' 2Nsdz 2 



d 2 d 

zP(z, t) — cv — (1 — z/z)zP(z, t) 



(18) 



where we have denoted xo by z for simplicity. The fact that the fittest class is lost whenever its size hits corresponds 
to an absorbing boundary condition for P{z, r) at z = 0. For such a one dimensional diffusion problem, the mean 



first passage time can be computed in closed form ( Gardiner 2004 1 and this formula has been used in ( GORDO and 
CharlesworthI 120001 IStephan et al] 119931) to estimate the rate of the ratchet. An accurate analytic approximation 
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to that formula has been presented by Jain (2008). For completeness, we will present an alternative derivation of 
these results that will help to interpret the more general results presented below. In the limit of interest, Nsz ^> 1, 
clicks of the ratchet occur on much longer time scales than the local equilibration of z. We can therefore approximate 
the distribution as P(z, r) « e~ 1T p(z) where 7 is the rate of the ratchet. In this factorization, p(z) is the quasi-steady 
distribution shown in Fig. [3]^, while 7 is the small rate at which P(z,t) looses mass due to events like the one shown 
in Fig. ^p- Inserting this ansatz and integrating Eq. (18) from z to 00, we obtain 



Id 

jP(X > z) = — — -w-zp(z) - az(l - z z)p(z) 
2Ns oz 



(19) 



where P(X > z) — dz'p(z'), which is ps 1 for z < z and rapidly falls to for z > z. To obtain the rate 7, we will 
solve this equation in a regime of small z<z, where the term on the left is important but constant, and in a regime 
z ^> (Ns)^ 1 , where the term on the left can be neglected. For the general discussion below, it will be useful to solve 
this equation for a general diffusion constant D(z) (here equal to z/2Ns), force field A{z) (here equal to az(z/z— 1)), 
and a constant C 



C 







D(z)p(z) + A(z)p(z) 



(20) 



with solution 



p{z) 



D(z) 



- P dv A(y) 

_g Jo "»D(i) 



P - C / dye 



rv dv i My') 



(21) 



Note that this solution is inversely proportional to the diffusion constant, while the dependence on selection is 
accounted for by the exponential factors. For z <C z, C = 7 and 2aNs J* dy (1 — z/z) ps 2aNsz, such that 



JINs 



p(z) PS 7- 



1 



z«z 



Q Z 



(22) 



where f3 is fixed by the boundary condition that p(z) is finite at z = 0. Note that 7 = p(0)/2Ns relates the rate of 
extinction to p(0). To determine the latter we need to match the z <C z regime to the bulk of the distribution z ps z. 
As can be seen from Eq. (22), the constant term C is unimportant in this regime (e 2JVsQZ ^> 1). Setting C = in 
Eq. (pIL we find 



p(z) 



liVs 



zNsa > 1 



(23) 



The integration constant (3 in Eq. (21 ) is fixed by the normalization. Since p(z) is concentrated around z = z and has 
a Gaussian shape around z, the normalization factor is simply 1/V27r<7 2 , where a 2 = z/2aNs is the variance of the 
Gaussian. The factor z/z corresponds to the l/D(z) term, scaled so that it equals 1 in the vicinity of z. Note that 
we have already calculated the variance of p(z) earlier, Eq. (13), and that consistency with this result would require 
that a is determined by Eq. (13). 

The two approximate solutions Eq. (22) and Eq. (23) are both accurate in the intermediate regime (Nsa)^ 1 <C 
z <C z, which allows us to determine the rate 7 in Eq. (22) by matching the two solutions. This matching implies that 
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\fzW& 



iNs 



(24) 



which agrees with results obtained previously (Jain 
and Ns of the rescaled model 



2008) 



Note that this rate only depends on the parameters A 
Since rates have units of inverse time, this expression has to be multiplied by s to 



obtain the rate in units of inverse generations. 

However, Eq. ( 24 ) does not describe the rate accurately, as is obvious from the comparison with simula tion results 
shown in Fig. Uk. The plot shows the rescaled ratchet rate 7 x y/n/ V zNsa 3 , which according to Eq. (24) should 
be simply exp"aiVse _A ], indicated by the black line. The plot shows clearly that the simulation results often differ 
from the prediction of Eq. (24) by a large factor. It seems as if a needs to depend on A, as we already noticed above 
when comparing the variance of p(z) to Eq. (13). In fact, fixing a via Eq. (13) improves the agreement substantially, 



but still does not describe the simulations quantitatively. 

The reason for the discrepancy is the time delay between 5k and z which we quantified by calculating the correlation 
between dk(r) and z(t + At). Hence we cannot use an approximation where Sk depends on the instantaneous value of 
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FIG. 4 The ratchet rate from simulation vs prediction. Both panels show the ratchet rate 7, rescaled with a prefactor to 
isolate the exponential dependence predicted by analytic approximations; A is color-coded. Panel A compares the simulation 
results to the prediction of Eq. (24 1, which is shown as a straight line. The approximation works only for a particular value of 
A, for otherwise the exponential dependence on Nxqs is not predicted correctly. Panel B compares simulation results to the 
prediction of Eq. (32 1, again indicated by the straight line. The exponential dependence of rate on iVsS^(O) is well confirmed 
by simulation results. 



z, but have to calculate 6k from the past trajectory of z. If the fittest class is the only one that is strongly stochastic, 
we can calulate dk(r) for a given trajectory z{p), p < t by integrating the deterministic evolution equations for Xk 
with k > 1 with z(p) as an external forcing. 

Eq. (17) now not only depends on z(t), but on all z(p) with p < r and cannot be mapped to a diffusion equation. 
Nevertheless, it corresponds to a well defined stochastic integral, known as a path integral in physics ( |Feynman and| 



Hibbs 19651, which is amenable to systematic numerical approximation. To introduce path integrals, it is useful to 



discretize Eq. (171 in time and express z{pi) in terms of the state at time = Pi — At and the earlier time points. 
For simplicity, we will use the notation Zi for z(pi). 



-1 = ArSk^xZi-x + 



-iAr 



Ns 



Vi-i 



(25) 



where 8ki-\ depends on all previous time points pj with j < i. In the limit At — > 0, this difference equation converges 
against Eq. (17) interpreted in the Ito sense since the z dependent prefactor of the noise term is evaluated at pi_i 
rather than at an intermediate time point between pi-\ and pi. We can express this transition probability P t (z\zq) 
between the initial state zq and the final state z = z m as a series of integrals over all intermediate states z% for 
< 1 < TO. 



~ m— 1 

P T (z\z ) = / dz i P/! i . T (z\{z j } j<m )PAT(zm-i\{z j } j<m - 1 ) ■ ■ -Pat(zi\zo) 



(26) 



Each of these infinitesimal transitions correspond to solutions of Eq. ( 25 ) with 7^ drawn from a standard Gaussian 
(ILau and LubenskyI 120071). Hence 



PAT(Zi\{Zj}j<i 



'Ns 



■ exp 



-Ns 



(&i - Zi-i - Arzj_i(5fci_i) 
2Atz,_i 



^/27rAtzi_i 

In the limit of many intermediate steps and small At, the transition probability can therefore be written as 

121 



/ ' Vz p exp 


-Ns j 




Jo 



dp- 



Zp5k pi 



2z„ 



= I Vz p e~ NsS ^}) 



(27) 



(28) 
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where T>z p is the limit of Y\T= 1 dzi(2nATZi-i/Ns)~ 1 / 2 known as the path integral measure, and we have replaced 
the discrete time index by its continuous analog. The path integral extends over all continuous path connecting the 
endpoints z and z T = z. The functional S\{{z p }) in the exponent closely corresponds to the "action" in physics 



(Feynman and HlBBS. 19651 which is minimized by classical dynamics. Here minimization of the "action" defines 
the most likely trajectory. Note that S\({z p }) depends on the entire path {z p } with < p < t, while the functional 
itself only depends on A. The strength of genetic drift appears as a prefactor of S\({z p }) in the exponent. 

The most likely path z* connecting the end-points points z and z T in time r can be determined either by solving 
the Euler-Lagrange equations or by numerical minimization, see below. Along with the functional, z* depends only 
on A. Given this extremal path, we can parameterize every other path connecting zq and z T as z p = z* + 5z p , where 
dzp vanishes at both endpoints (Szq = Sz T = 0). Denoting the minimal action associated with z* by S^(z T , zq), we 
have 

P T (Z T \Z ) = e -NsSl(z T ,z ) J e -Ns SS \({8z p } ,z T ,z ) = R - N *S\ (z T ,*„) ( 2 9) 

where Af~ l factor is equal to the integral over the fluctuations, which in general depends on z*. The prefactor 
Ns in e ~ Ns ss ({ Sz p}> z t,zo) i m pli e g that deviations from the optimal path are suppressed in large populations. If 
6S({5z p }, z T , Zo) is independent of the final point z T , Af can be determined by the normalizing P t (z t \zq) with respect 
to z T . In the general case, calculating the fluctuation integral is difficult, and we will determine it here by analogy to 



the history independent solution presented above Eq. (18|-(24). 

If the stochastic dynamics admits an (approximately) stationary distribution, P T (z T \z ) becomes independent of t 



and zq and coincides with the steady state probability distribution p(z). It therefore becomes the analog of Eq. (23 1, 
which for arbitrary diffusion equations is given by the inverse diffusion constant (the prefactor multiplied by an 

exponential quantifying the trade-off between deterministic and stochastic forces. In this path integral representation, 
the exponential part is played by e~ NsS ^ z \ where S^(z) is a function of the final point z and A only. The prefactor is 
independent of the selection term and can hence be determined through the analogy to the Markovian case discussed 
above 

p(z) a - . 1 e - yjg *W Nsz > 1 (30) 
where the normalization is obtained by assu ming an approximately Gaussian distribution around the steady state 



value z and the variance a 2 is given by Eq. (13). Note that this solution is not valid very close to the absorbing 
boundary since this boundary is not accounted for by the path integral, at least not without some special care. As 
in the history independent case discussed above, this approximate distribution should be thought of as the time 
independent "bulk" distribution in P(z,t) = e ,T p(z). To determine the rate extinction rate 7, we again need to 
understand how probable it is that a trajectory actually hits z — 0, given that it has come pretty close. 

To this end, we need a local solution of Eq. (Tl7| in the boundary layer z <§; z as already obtained for the history 



independent case in Eq. (22). Once a trajectory comes close to z = 0, it's fate, i.e., whether it goes extinct or returns 
to z ~ z, is decided quickly. Hence we can make an instantaneous approximation for Sk which does depend on the 
past trajectory, but for the time window under consideration it is simply a constant, a, yet to be determined. Having 



reduced the problem to Eq. (22) we can determine a, and hence 7, by matching of the boundary solution to Eq. (30) 
in the regime (Ns)^ 1 « z « z where both are accurate. The matching condition is 

Z 1 e -N s (S*(.0)+zd z S*(z)\ z=o ) _ J_ e 2aNsz /g^ 

z V2ira 2 az 
which determines a and 7 by the matching requirement 2a — — d z S^(z)\ z=0 and 

\d z Sl(z)\ze- N ^ ^ S* x (0) c _ NsSU0) 
2V2^ ~ V^a 2 



where we approximated \d z S^(z)\ z= o m S^(0)/z. The variance a 2 is given by Eq. (13) and depends on A and Ns 
as a 2 — zC t 2 (X)/Ns. Note that 7 is in units of s and needs to be multiplied by s for conversion to units of inverse 
generations. In contrast to Markovian case above, the variance of the "bulk" is no longer simply related to strength 
of selection near the z = "boundary" . 

Since we don't know how to calculate S^(0) or the most likely path z* analytically, we determined discrete ap- 
proximations to z* numerically as described in section Model and Methods. Examples of numerically determined 
most likely path and the corresponding trajectory of the mean fitness are shown in Fig. [5] for different values of A. 
Generically, we find a rapid reduction of the size z of the size of the fittest class such that the mean fitness has only 
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FIG. 5 Panel A: The most likely path Xq(t) to extinction of the fittest class and the concomitant reduction of the mean fitness 
for different A are plotted against time. Times are shifted such that Xq(0) = The inset shows the mean fitness Sk(r) plotted 
against Xq(t) for different values of A. Panel B: "Haigh's factor" q(A) = S^(0)e A as a function of A determined numerically. 



partially responded. The inset of Fig. [5j\ shows how changes in mean fitness 8k are related to z for different A. For 
large A, the mean fitness changes only very slowly with z, which increases the probability of large excursions and 
hence the rate of the ratchet. 



This numerically determined minimal action S* x (0) together with the approximation Eq. ( 32 1 describes the rate 
of the ratchet, as determined in simulations, extremely well. Fig. [4j3 shows the same simulation data as Fig. [3]^, 
but this time rescaled by yj^a 2 (A) / <S^(0) as a function of NsS^(0). After this rescaling, we expect all data points 
to lie on the same curve given by exp[— 7VsS^(0)], as is indeed found for many different values of it, s and N with 
A = u/s ranging from 1 to 30. Note that the vertical shift of the black line relative to the data points depends on the 
prefactor, which we have approximated. Hence we should not expect agreement better than to about a factor of 2. 
The important point is that the exponential dependence of the rate on NsS^Jfi) is correctly captur ed by Eq. (|3 2|). 



Previous studies of Muller's ratchet suggested that the rate depends exponentially on aNse (Jain, 2008). To 



relate this to our results, we determined "Haigh's factor" a numerically from 5^(0) and plotted it in Fig. |5J3. We 
find that a(A) drops from around 0.8 to 0.3 as A increases from 1 to 30. The previously used values 0.5 — 0.6 for a 
correspond to A ~ 6. Using cv(A) = S^(0)e A as shown in Fig. [5^3, we can recast Eq. (32) into its traditional form and 
undo the scaling with s. In units of generations, the mean time between clicks is given by 



T, 



click 



2.5C(A) 

_ Nsa(X)e 



(33) 



where £(A) is determined by Eq. (13) and the factor 2.5 is introduced to approximate the part of the prefactor that is 
independent of N, s or A. A direct comparison of this expression with simulation results is shown in supplementary 
figure 1. 



III. CONCLUSION. 

The main difficulty impeding better understanding of even simple models of evolution is the fact that rare events 
involving a few or even single individuals determine the fate of the entire population. The important individuals are 
those in the high fitness tail of the distribution. Fluctuations in the high fitness tail propagate towards more mediocre 
individuals which dominate a typical population sample. 

We have analyzed the magnitude, decay, and propagation of fluctuations of the fitness distribution in a simple 
model of the balance between deleterious mutations and selection. In this model, individuals in the fittest class evolve 
approximately neutrally Fluctuations in the size of this class propagate to the mean, which in turn generates a delayed 
restoring force opposing the fluctuation. We have shown that the variance of the fluctuations in the population n of 
the top bin is proportional to uq/s and increases as log A with the ratio A of the mutation rate u and the mutational 
effect s. Fluctuations of no perturb the mean after a time ~ s _1 log A. These two observations have a straightforward 
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connection: sampling fluctuations can accumulate without a restoring force for a time s 1 log A. During this time, the 
typical perturbation of the top bin by drift is ~ \J uqs^ 1 log A and hence the variance is sa uqs' 1 log A. We have used 
these insights into the coupling between no, the mean fitness, and the resulting delayed restoring force on fluctuations 
of n to approximate the rate of Muller's ratchet. 

The history dependence of the restoring force has not been accounted for in previous analysis of the rate of Muller's 



ratchet (Gordo and Charlesworth 2000 Haigh 1978 Jain 2008 Stephan et al. 1993) who introduced a 



constant factor to parameterize the effective strength of the selection opposing fluctuations in the top bin, or |WAXMAN 



and Loewe (2010), who replaced all mutant classes by one effective class and thereby mapped the problem to the 



fixation of a deleterious allele. We have shown that to achieve agreement between theory and numerical simulation one 
must account for the delayed nature of selection acting on fluctuations. Comparing our final expression for the ratchet 



rate with that given previously (GORDO and CHARLESWORTH 2000 Jain 2008 Stephan et al. 1993), the history 
dependence manifests itself as a decreasing effective strength of selection with increasing A = u/s. This decrease is 
due to a larger temporal delay of the response of the mean fitness to fluctuations of the least loaded class. History 
dependence is a general consequence of projecting a multi-dimensional stochastic dynamics onto a lower dimensional 
space (here, the size no of the fittest class). Such memory effects can be accounted for by the path- integral formulation 
of stochastic processes which we used to approximate the rate of Muller's ratchet. 

Even though the model is extremely simplistic and the sensitive dependence of the ratchet rate on poorly known 
parameters such as the effect size of mutations, population size, and mutation rate, precludes quantitative comparison 
with the real world, we believe that some general lessons can be learned from our analysis. The propagation of 
fluctuations from the fittest to less fit individuals is expected to be a generic feature of many models and natural 
populations. In particular, very similar phenomena arise in the dynamics of adapting populations driven by the 
accumulation of beneficial mutations flCOHEN et al. 2005 Desai and Fisher 2007 Hallatschek 2011 Neher 
et al.\ |2010| |Rouzine et al.\ |2008| |2003| |TsiMRlNG et a/^ |1996[ ). The speed of these traveling waves is typically 



determined by stochastic effects at the high fitness edges. We expect that the fluctuations of the speed of adaptation 
can be understood and quantified with the concepts and tools that we introduced above. 

Populations spread out in fitness have rather different coalescence properties than neutral populations, which are 
described by Kingman's coalescent ( Kingman 1982[ ) . These differences go beyond the familiar reduction in effective 



population size and distortions of genealogies due to background selection ( Charlesworth et al. 1993 



Woodcock 



1995 Walczak et al. 2011). The most recent common ancestor of such populations most li 



HiGGS and 



lely derives 

from this high fitness tail and fluctuations of this tail determine the rate at which lineages merge and thereby the 
genetic diversity of the population (|Brunet et al] |2007| |Neher and Shraiman| |2011| |Rouzine and Coffin] |2007 1 . 



Thus, quantitative understanding of fluctuations of fitness distributions is also essential for understanding non-neutral 
coalescent processes. 

Generalizing the analysis of fluctuations of fitness distributions to adapting "traveling waves" and the study of their 
implications for the coalescent properties of the population are interesting avenues for future research. 
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Appendix A: Supplement 

To analyze the covariances of different fitness classes, we begin with Eq. (7) of the main text, which expresses 8xk{j) 



in terms of eigenvectors 8xk{T) = y~) . ip^ Oj(t). Projecting on the left eigenvectors then results in equations for Oj(t) 



± aj (r) = -j a, (r)+5>?V|£»&(r) ( Al) 

k 



where ^(t) are uncorrelated Gaussian white noise terms with {rjk(j)rji{T')) = 5kiS(r — r ). Since each noise term r]k 
contributes to all a_,-, the noise induces correlated fluctuations of the dj{r), which we need to understand in order to 



analyze the fluctuations of the fitness distributions. The inhomogeneous Eq. ( Al ) has the solution 



,.'.'■'./■'■'■■ 

k ' 

The autocorrelation function of the loadings of different eigendirections separated by At in time is therefore given by 

dT /g-i(T-T')-j(T+AT- T ') ^ 9k 9k k (A3) 



Ns 

k 

k ,( Pk' x k 



e -jAr rh ^) rh U) 

~ i+ j ^ Ns 

J k 

where we have used (r)k{T)r)i(T')} = SkiS(r — t'). 

Correlation functions no and the mean 

To calculate the variances and covariance of xq and the mean fitness, we express them in terms of the eigenmodcs 
aj (r) (i>0) 

5x (t) = ^V^ ) «,(r) = - e - A ^a J (r) (A4) 

j>0 j>0 

<5fc(r) = - ^fc a j( T ) = - E K^k-j -Xk)aj(r) = -^jai(r) (A5) 

j>0,k j>0,k j>0 

The auto-correlation of xo 

The auto-correlation of xq is given by 

p -jAt A^Art*, 
e 9k Vk x k 

i + j 2— 1 Ns 

J k 



e- 2X 


E 




i,j>0 






Ns 


E 


i,j>0 


e- A 




Ns 





l+ i i~o U - - k)\kl 



(j-k)Hi-k)\k\ 

i,j>0 fc=0 w ' y ' 
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Let us focus on the triple sum inside the integral and simplify it by introducing a = — \e z and b = — Ae 2 At . 
Furthermore, let us look at the i = j and the i 7^ j contributions separately. The diagonal contribution (i = j) is 



V 'ft'V A ^ = V V (*by- k (ab) k \- k 

2-^ 2- (i - k)\(i - k)\k) 2^2^ k\((i~k)\) 2 

i=0 fe=0 v /v ' i>0 fe=o vv '' 



EE 

fc>0 

E 



{aby- k {ab) k \- k 



V W 

m-kW "^(i!) 2 



fc>0 i>k KX ' ' i>Q 

{ab) k \- k ^ (ab) n 



(A7) 



E , , xo + Jo(—L2Vab) — 1 using n = i — k 
IT?.' 



fc! ^„ (n!) 2 



fc>0 ri>0 

= (e ab/A - l)J (-i2\/a6) + J (-i.2\/a6) - 1 
= e ab / A J (-t2Val^) - 1 

where J n (z) is the nth Bessel function of first kind, and l = \] — 1. When evaluating the off-diagonal contribution, we 
will encounter terms like 



(A8) 



(ab) k _ (ab) k 1 _ J m {2tVab) 1 

^^fc!(fc + m)! k\(k + to)! ml (i\/ab) m m\ 



The off-diagonal contribution can be further split into the parts i > j and i < j which can be evaluated as follows: 
^ 2^ kUj - k)\(i -k)\ 2^Z^ UU 2-, kl(i + (7 - i) - k)\(i - k)l 

0<i<j k>0 w ' y ' i>0 j>i k>0 v w ' ' v ; 



i>0 m>0 

EEE 



A" 



_ o fc!(i + to - fc)!(i - fc)! 



using m = j — 1 



fc!(i + m - k)\(i - k)\ 2—/ 2—, u + m )n\ 

k>0m>0i>k y I \ 1 m>0 i>0 v ' 

n+kum+n+k\—k i n h\i 

y^ y^ y^ a g ^_ , y^ fe ™ y^ l a "J 

M(n + m)\n\ 2—, 2^fi + m \\n 

fc>0m>0n>0 V 7 m>0 i>0 v 1 



using n = i — k 



(A9) 



a k b m+k X -k j m (2,Vab) y- & „ 
2—i 2—, h\ ( L Jab\^ 2—, 

fc>0m>0 \L\UU) m>Q 



J m (2tva&) 1 
(t\/afe) m m! 



^-v ^ a k b k b m / 2 a- m /' 2 \- k J ro (2tVa6) ^ 



m>0 fc>0 



m>0 



fc! 

l/2 



(0' 



(-l)™/2 e +i 

m>0 v 7 v 7 



(-0 m ^m(2^Vafe) - e" + 1 



The off-diagonal terms for i > j is obtained by interchanging a and 6 such that the full off-diagonal contribution is 

m / 2 ,„.m/2 



a ab/X 



E 

m>0 



(~i) m J m C2tVab) -e a ~e b + 2 



(A10) 



Next, we use the definition of the generating function of the Bessel functions (Gradshteyn and Ryzhik (20071, 
8.511) 



e Ut-t-^ = J o{z) +Y / (t m + (-t)- m )Jm(z 



m>0 



which turns the off-diagonal contribution into 

s a6A ( e ^~ b (^ + ^) - J (2iVab) )~e a ~e b + 2 



(All) 



(A12) 
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Combining the diagonal and off-diagonal contributions and substituting a and b, we find for the integrand in Eq. ( A6 ) 



.ab/X+Vab^+Jf ) _ e a_ e b + 1 = e Ae — -Ae^-Ae— _ e -A e - _ fi -Ae— + j (A1 gj 



The auto-correlation of xo is therefore given by 

— A poo 

(x (r)x (r + At)) = [ dz (f*-—*\-K*-'+*—*) e -*- - e -A.— A * + X 



o 



— ^ e A(?2e ~ AT - Ae ( 1 + e ~ A ' r ) _ e - Ae _ e - ASe_AT + 



(A14) 



Auto-correlation of the mean fitness 

The autocorrelation function of the mean is defined as 

(Sk(T)Sk(T + At)) = ^ ij(oi(r)oj(r + Ar)) 



»,J>0 

min 



1 / j_ 7' -2(J+i) \i+7' -;'Ar \ "* V 1 i A 



(j -k)Ui- k)\k\ 

i,j>0 k=0 w ' v ; 

where the last line is to be evaluated at v = ji = 1. Hence the problem is reduced to the one already solved with 
a = — /jAe~ z and b = — v\e~ z ~ Ar . We find 

(5k{r)5k{r + Ar)) = ^ t dee^e™ 2 ^-^'^ 6 (6 + \6 (Oe^ - l) (9 - 1)) (A16) 
Ns J 

Cross-correlation of x and the mean fitness 

When calculating the cross-correlation between xq and the mean fitness we have to distinguish the cases where 
xq precedes the mean fitness and vice-versa. Otherwise, the calculation proceeds almost unchanged from the cases 
discussed above. 



[Xq 



( T )Sk(T + Ar)) = J2 3{ai{T)aj{T + At)) 



i,j>0 

1 roc min(ij) (_ 1 ^-fc ( A17 ) 

with a = -\er z , b = -v\e~ z - Ar if Ar > and a = -Ae~ 2+Ar , b = -v\e~ z if Ar < 0. The result is 

A fe- A - tidefte-l^-^^-^+e-^ + e-^™) Ar>0 

(fa0(T)ft(T + AT)) = -^ 1 / V At 2 / At \ ' (A18) 

V W 1 " Ns U^d6((e Ar 6-l)e e Txe - x ^+ e ) e + e- xe ) Ar < 
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FIG. 6 The approximation of the mean time between clicks of the ratchet is accurate over a large range of parameters if 
Nsa(\)e~ x is large compared to one. Nsa(\)e~ > ' determines whether the clicks of the ratchet are far apart compared to the 
relaxation time of the distribution and is indicated as the color of the data points. The condition jVsa(A)e _A > 1 is violated 
for the fastest clicks shown, resulting in the deviation of the dark blue points. 



